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I. INTRODUCTION 



Inversion of integral transforms which values are known only approximately is a long- 
standing issue. The situation is unfavorable when the solution that corresponds to an 
approximate input differs from the true solution in a way that it includes extra narrow 
peaks or quick oscillations. In these cases small changes in the input correspond to changes 
in the solution that are not small. 

But in such cases the numbers of maxima or minima of the approximate solution and of 
the true solution would differ from each other. Therefore, to cure the situation it is suggested 
in the present paper to seek for the approximate solution in the class of functions having 
the same number of maxima and minima as the true solution has. In general this number is 
not known but often it is not hard to guess it as discussed below. A general ansatz for the 
solution compatible with a prescribed number of maxima and minima is proposed below for 
performing inversion. 

In Sec. 2 the proposed inversion method is described. It is applicable to a wide class of 
integral equations of the first kind. 

In particular, the techniques are aimed to increase ability of the integral transform ap- 
proach used to calculate amplitudes and inclusive spectra of perturbation-induced reactions 
in the framework of nuclear physics. An outline of that approach is presented in Sec. 3. 

In Sec. 4 the techniques are tested with exactly solvable examples. The Laplace, Stieltjes 
and so called Lorentz [T] integral transforms are inverted. 

In Sec. 5 the issue of convergence of the method to the true solution is investigated. The 
important point is that the solution is sought for in a restricted class of functions. This 
ensures the convergence. (A different case with such features exists, and in this case the 
solution belongs to a known compactum, see e.g. [2].) 

II. INVERSION METHOD 

The integral equation 



is studied. All the quantities are real, the solution / is smooth and unique. The K operator 
is assumed to be continuous at the norm definitions specified below. The exact $ is not 





(1) 
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known so that one is forced to deal with an approximation to $ denoted below $appr- In 
the case of few-body calculations performed in the framework of the approach outlined in 
Sec. 3 the accuracy of an input $appr is normally at a per cent level. (The notation for the 
lower integration limit in Eq. ([T]) is used in connection with the problems addressed in Sec. 
3 and 4. In the same connection the upper integration limit is set to be infinite.) 

As it is known the solution / sought for may be unstable with respect to small changes of 
the input. In the literature, various regularization procedures were put forward to suppress 
the instability. (A rather complete review can be found in Refs. ^ E], see also e.g. [HE].) 
In particular, within the approach outlined in Sec. 3 the following regularization procedure 
was always used so far. The solution was approximated by the expression 

N 

fN{E) = J2cnME,a) (2) 

n=l 

where {v?n} is a set of basis functions. The linear parameters Cn and the non-linear parameter 
a were found via fitting the quantity Kf^ to $appr- The regularization was realized with 
the choice of the N value. This value should be not too high to exclude unstable behavior. 
At the same time it should be not too low so that the approximation of the true f{E) with 
the Eq. g type expression would be sufficiently accurate. (See e.g. Ref. p| for further 
details.) The problem is that the two conditions are compatible with each other only to a 
certain degree. If the accuracy in an input $ is not high enough reasonable inversion may 
become impossible or tricky. The issue of the proper choice of the regularization parameter 
does not arise in the method of the present paper since no regularization is required here. 

The techniques proposed below are designed to increase the accuracy of inversion at a 
given accuracy in $. It should be noted in this connection that when Eq. (|2| is used another 
source of instability exists besides uncertainties in the input. This instability is related to 
arising ill-conditioned systems of linear equations for the expansion coefficients from Eq. (|2]). 
The influence of corresponding round-off errors on the solution becomes considerable when 

increases even if the exact input is employed. In the framework of the method described 
below the round-off errors are probably less important. 

It is shown in the last section that when small uncertainties in an input are of non- 
random nature and the round-off errors in calculation are not substantial no other sizable 
deviations from the true solution besides narrow peaks of small strength may exist in the 
problem. Adopting that this is the case let us impose the requirement that the solution 
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sought for includes a given number of maxima and minima. Then one may conclude, see 
also the last section, that if the guess as to their number is correct then those narrow peaks 
are excluded and correspondingly the approximate solution is close to the exact solution 
"almost everywhere". 

The only possible exception are the points of maxima and minima of the solution where 
extra narrow peaks might appear. If this improbable situation takes place these peaks are 
in general to be simply removed. Indeed, true narrow peaks may appear only when there 
exist specific physical reasons for this, such as resonant states. When the above requirement 
on the number of maxima and minima is imposed no regularization is required. Below the 
solution is sought for in the class of functions satisfying this requirement. 

In many cases incorrect guesses as to the number of maxima and minima of the solution 
are easily rejected. Some maxima or minima may be very weakly pronounced. In such 
a case it may occur that approximate solutions with different numbers of maxima and 
minima are very close to each other and lead to the fits of a similar quality. Each of these 
approximate solutions is acceptable. Apart from this case, if the inversion is performed 
under the condition that the number of maxima and minima is less than the true one this 
necessarily should lead to a fit of a lower quality. And if the inversion is performed under 
the condition that their number is higher than the true one this should lead either to a fit 
of lower quality or to unrealistically narrow extra peaks in the solution with peaks positions 
and amplitudes being unstable. 

The inversion procedure is as follows. In the problems addressed in Sec. 3 the threshold 
behavior of f{E) at E ^ Ethr is known. ^ In particular, one has f{Ethr) = 0. Also f{oo) = 0. 
Let us rewrite Eq. ([T]) in terms of df /dE = f'{E) integrating by part: 

/ k{a,E)nE)dE = ^{a). (3) 

Ethr 

{K{a,E) = — J K{a, E)dE.) Eq. (jsj) supplemented with the conditions 

f{E,f,r) = 0, /(oo) = (4) 
is equivalent to Eq. ([T| at inputs $ that lead to the values Q. To satisfy the conditions Q 



^ This behavior is deduced from the known behavior at E.y — Ethr of the matrix elements entering Eq. ( 17). 
In this hmit they include the E-y dependence as a factor. 
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let us impose the requirement 

POO 

/ f'iE)dE = /(oo) - /(Ei,,) = (5) 

and subsequently find the solution to Eq. (1) e.g. as f{E) = f'{E')dE'. The solution 
to Eq. ([3]) with the condition ^ is unique. Indeed, otherwise there would exist more than 
one solution to Eq. ([T]) constructed in the above way. 

To be specific let us consider below the most frequent case when at high E values the 
solution f{E) decreases as power of E. (In the case of problems addressed in Sec. 3, one 
can see that in order to lead not to a power but to an exponential decrease of f{E) the 
inter-particle interaction in the coordinate or momentum representation should be analytic 
in all the coordinates, or momenta, respectively. This is not the usual case.) Let f{E) 
have precisely N maxima and minima. Let (thr{E) be a monotonous factor reproducing the 
threshold behavior of f{E) at — > Ef^r, behaving as power of E (e.g. as a constant) at 
E —¥ 00, and arbitrary otherwise. The following ansatz for f'{E) is suggested. 



fiE) = CCuAE) 



N 

(6) 



[{E/E) + 1] 



j=i 

Here CthriE) = d(thr/dE, E and /3 are parameters, 7(-E) is a smooth function finite both at 
E = Ethr and at E ^ 00, and e.g. 'j{Ethr) = 0. The parameter C determines the overall 
normalization. This expression is the most general one for the derivative of a function that 
has precisely maxima and minima, a given threshold behavior, and a power decrease at 
-E — > 00. (Infiection points do not require a special consideration. A given power decrease 
is provided by the choice of the parameter /3.) The function 'j{E) may be taken e.g. in the 
form 

AE/E ^ c, 
' {AE/E) + lf^^[{AE/E) + l]^ 

with AE = E — Ethr- Eq. ([T]) allows the description of next-to-leading terms in f'{E) when 

both E — > Ethr and E — t- 00. The fitting parameters are then C, {-Ej}, E, /3, and {cjt}. 

To satisfy the requirement ([s]) it is convenient to express one of the Ei values in ([6]) in 

terms of the other fitting parameters. In certain cases also the sum-rule condition 



/•oo /»oo 

/ Ef{E)dE= / f{E)dE = S 



(8) 



may be imposed where the value of S is known, see Sec. 3. From this condition it is 
convenient to express C in terms of the other fitting parameters. 



The remaining parameters are to be determined from the least-square fit procedure 

1 1 ^appr - Kflj 11 = 11 ^appr - K/m \ \ = Hlin. (9) 

Here fj^^ is an approximation to /' of the form (j6]) with M parameters retained, /m is the 
corresponding approximation to /, and the norm is defined as 

1/2 



w{a)F\a)da 



w{a) > 0. (10) 



(When the equation of Eq. ([T]) type is considered with a and b being the lower and the 
upper integration limits and when /(a) ^ and/or f{b) ^ one may proceed in a similar 
way. One may rewrite the equation in the form of the equation for the derivative f'{E) 
and from that equation one obtains f'{E) = p{E)f{a) + q{E)f{b) + r{E) where p, q, and 
r are known functions. Then it is easy to get the relation f{E) = A{E) + \B{E) where 
A may be chosen to be /(a), /(&), or /(c), c being an arbitrary point, and A and B are 
the corresponding known functions. If there is no point c at which /(c) is known then the 
parameter A is to be determined from the initial integral equation.) 

In simple cases the main features of the solution f{E) are determined by its threshold 
behavior, its decrease at large E, and by positions and amplitudes of its maxima and minima. 
For example, let us consider here the one-maximum or one-minimum case. The mentioned 
features of the solution may be reproduced in the framework of the ansatz ([6]) without 
the factor exp[7(i?)]. Therefore one may think that the solution is rather insensitive to the 
values of the parameters determining 7(-E). In addition, in many cases the dependence of the 
quantity || $appr ~ -^/mII t)e minimized on the parameters is rather smooth. Then a good 
strategy to find the global minimum of Eq. (|9| would be in the first stage to disregard the 
factor exp[7(ii^)] in Eq. ^ and to seek for the minimum with respect to other parameters 
on a grid that includes sufficiently large number of points. Grids where one grid is put 
inside another one may be used. The values of the parameters thus obtained may be used 
as starting values at a subsequent conventional minimization. 

To carry out the latter minimization a good choice are special codes for the least-square 
fit with non-linear parameters (e.g. Ref. [7|, Sec 15.5). Codes that use derivatives are 
preferable since the derivatives of e.g. the expressions ^ and ([T]) with respect to fitting 
parameters are readily calculated. 
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At calculating the norm in ([9]) one replaces the integral with an integral sum. A good 
approximation to the solution may be obtained even if the number of terms in this sum is 
insufficient to reproduce the integral. 

Practical criteria of the quality of an approximate solution Jm thus obtained are both 
its stability with respect to increase of the number M of fitting parameters at a given $appr 
and its stability with respect to variations of ^appr- If a realistic estimate of 11$ — $appr || is 
available, say ||$ — $appr|| > ^o, then the criterion of stability with respect to increase of the 
number M of fitting parameters may be replaced with the condition that the quality of the 
fit, i.e. the value of \\KfM — ^appAl, is comparable with eo. 

After this work has been completed I got to know about the work [8] where an approach 
having some common features with the present one was considered. The techniques of 
Ref. [8] are different from those of the present paper. The authors consider segments of 
monotonicity of the solution i.e. segments of a constant sign of its derivative (and also of 
constant curvature of the solution). Discretizing the problem they reduce the corresponding 
minimum condition to the problem of the quadratic programming with linear constraints. 
However, no systematic way was given in Ref. [8] to find the borders of the monotonicity 
segments i.e. optimal positions of maxima and minima. Thus, unlike the present method, 
in Ref. |8j the whole problem seems to have not been solved. Besides, the issues related to 
the appearance of unrealistically narrow peaks discussed above (in relation to a guess as to 
the true number of maxima and minima) and in the last section (in relation to convergence 
of the method) are not considered in Ref. |B]. 

It may also be noted that the techniques of Ref. [8j require finding the minimum with 
respect to many variables (like = f{Ei)) with constraints. While in the present techniques 
finding the minimum with respect to only a few fitting parameters without constraints is 
required. Furthermore, the behavior of the solution when it approaches the integration limits 
is not reproduced exactly in the techniques of Ref. [S]. Exact in-advance-reproduction of 
this behavior in the present techniques increases an overall accuracy of the solution. 
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III. INTEGRAL TRANSFORM APPROACH FOR CALCULATING REACTIONS 
IN QUANTUM MECHANICS 



The material of the present section is used below to generate realistic inputs ^appr- Some 
features of the ansatz of the preceding section for the solution f{E) are also related to this 
material. 

The techniques of the preceding section are designed, in partcular, for applications in the 
framework of the approach that is reviewed below. Only a brief outline of the approach is 
contained here. More details can be found in the reviews [BJ [S] . The approach is advan- 
tageous, in particular, in problems with many open channels of various nature i.e. when 
energy is not low. 

Its main features are the following. The dynamics calculations to be performed are 
bound-state type calculations. In the course of calculations there is no need to consider 
reaction channels, as well as reaction thresholds. Reaction channels and thresholds come 
into play at merely the kinematics level only after a dynamics calculation is done. 

Continuum spectrum states never enter the game. In place of them, "response-like" 
quantities of the type 



R{E) = J2{Q'\'^n){^n\Q)S{E - E„) + ^ d^{Q'\^^){^^\Q)S{E - E, 



(11) 



are basic ingredients of the approach. Here are bound states and '^^ are continuum- 
spectrum states. They represent a complete set of eigenstates of the Hamiltonian H oi a. 
problem. The subscript 7 denotes collectively a set of continuous and discrete variables 
labeling the states which is symbolized in the summation over integration notation. The 
normalizations (\E'„|\E'„/) = 5n,n' and (\E'^|\E'y) = 5{pf — 7') are assumed so that 

Y,\^>n){^n\+^d^\^,){^>,\=I, (12) 

I being the identity operator. 



In the method discussed the quantities R{E) of Eq. (11) are obtained not in terms of 
the complicated states ^'-y entering their definition but via a bound-state type calculation. 
And reaction observables are expressed in terms of R{E) as quadratures. 

Let us first explain the latter of these points. Consider strong-interaction induced re- 
actions. Denote A(t)i{E) and A(t)f{E) the antisymmetrized "channel free-motion states". 
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Here the subscript i refers to the initial state of a reaction, the subscript / refers to final 
states of a reaction, (f)ij{E) are products of fragment bound states and of factors describing 
their free motion [TU] , and A denotes the operator realizing antisymmetrization with respect 
to identical particles Denote (f)i{E) = A{H - E)(j)^{E) and (f)f{E) = A{H - E)(pf{E). 
One has 0, = AV['^'^(l)i and 0/ = AVp''^(j)f, where V[j'^ are interactions between fragments 
in the initial and final states. Here it will be assumed that these interactions are of a short 
range so that the long-range inter-fragment Coulomb interaction is disregarded.^ The T 
matrix determining the reaction rates is [TU] 

T;, = Tj^- + (MmH -E- ie)-'\UE)), (13) 
e — )■ +0. Here Tj^"^"' is the simple Born contribution, 

= (0/10.) = (0/10.), 



and the main problem consists in calculating the second contribution in (13) that includes 



the Green function {H — E — ie) ^ . This contribution may be represented as 

j dE' Re{E'){E' -E- le)-^ (14) 



where 



Re{E') = J2{ME)\^n){'^n\4>i{EME' - E„) + ^d^{^f{E)\m,){m,\UEmE' - E,). 

(15) 



The quantity (15) is just of Eq. (11) structure (with the E ^ E' replacement). Thus, 



indeed, to calculate matrix elements of the T matrix it is sufficient to have quantities of this 



structure. Once they are available the integrations (14) are readily done. 

In order to calculate a perturbation-induced reaction amplitude (\l/j|0|\l/o) where O is 
a perturbation and \l/o is an unperturbed initial state the same is to be done with the 
0j — )■ 0\E'o replacement in the above relations. 



Now let us explain the above mentioned point on calculating the Eq. (11 ) type quantities. 
It should also be noted that such quantities may be of interest themselves representing 
observable response functions for inclusive perturbation-induced reactions. Let us rewrite 



This restriction can be weakened or removed. This is done in part in 
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Eq. (|TT} as 

R{E) = J2 RnS{E - En) + f{E), Rn = {Q'\^n) (^n|Q), (16) 

f{E) = ^ d^{Q'\^,){^,\Q)6{E - E,). (17) 
Calculation of the bound-state contributions Rn can be done separately, see also below. The 



contribution (17) includes an integral over few- or many-body continuum states that 
are very complicated except for low energies, and the problem just lies in calculating this 
contribution. If Ethr denotes the threshold value for continuum state energies then f{E) is 
different from zero at Ethr < E < oo. 

An easy task is the sum-rule calculation. Using Eq. ( [I2| ) one gets 

POO 

/ fiE)dE={Q'\Q)-J2Rn. (18) 

J Ethr n 

Obviously, the quantity (18) does not allow reconstruction of f{E) itself. To achieve this 
goal, let us consider "generalized sums" of the form 

POO 

/ K{a,E)f{E)dE. (19) 

Ethr 

Using Eq. (12) one obtains "continuous sum rules" 

K{a, E)f{E)dE = ^ dj{Q'\^,)K{a, E,){^,\Q) 

= {Q'\K{a, H)\Q) - J2 RnK{a, E^) (20) 



where as above H is the Hamiltonian of the problem and Rn are defined in Eq. (16). If one 



is able to calculate the quantity {Q'\K[a,H)\Q) entering Eq. (20) then one comes to the 
integral equation ([!]) for f{E) with 

$(a) = {Q'\K{a, H)\Q) - J2 RnK{a, En). (21) 

n 

And at proper choices of the kernel K one can completely reconstruct f{E) from this equa- 
tion. 

The presented approach to calculate reactions has been introduced in [12].^ The trans- 
forms with the kernels K{a, E) = {E — cr)~^ where p = 1 and 2 were employed. These are 



^ For observable responses R{E), i.e. in case of inclusive perturbation-induced reactions, a bound-state 
type calculation of their integral transforms has been suggested in in case of the Stieltjes transform 
and in |14j in case of the Laplace transform. Inversions of the transforms were not considered in those 
works. 
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the Stieltjes transform and the generahzed Stieltjes transform. Here a is chosen taking real 
values lower than the continuum spectrum threshold and ranging outside neighborhoods 



of the discrete spectrum values E^- In accordance with Eq. (21) the input in the integral 
equation is 

Denoting [H — cr)^^Q = \E' and {H — a)^^Q' = \E'' this can also be written in the two re- 
spective cases as 

= (g'l^) - V— ^ where (H - a)^ = Q, and (23) 

$(a) = (^'1^) 5^ where {H ~ a)^' = Q' , {H-a)^! = Q. (24) 

The states ^ and \1'' are localized. Therefore the inputs $(cr) are indeed calculable with 
bound-state type methods. 

The transform with the "Lorentz kernel" [T] was intensively used. The procedure of 
Eq. ^ was used for the inversion. The kernel can be written as 

K{a,E) = [{E-a*){E-a)Y' (25) 



where a = aR + iaj is now complex. The quantity $(cr) obtained in this case is of Eq. (24) 



form with the replacement {En — a) ^ — t- [{En — o'*){En — cr)] ^. One can also obtain the 



Lorentz input $(cr) from the dynamics equations that, like the Stieltjes case (23), involve 
only the initial-state source term Q. For this purpose one rewrites the above mentioned 
$(o") in the form [9] 

where ^^ = {H - a)-^Q, = {H - a*)-^Q. (26) 

Both ipi and 'i/'2 are calculated from the initial state of a reaction. Final states enter here 
via the known quantity Q' i.e. as quadratures. 



In the cases (23) and (24) it is convenient to calculate Rn as the limits of the expression 
— {En — (j){Q'\'^{a)) and of the expression — {En — cr)^(\E''(cr)|\^((j)) respectively at cr tending 
to En- Here and \I'' are the solutions to the corresponding inhomogeneous equations. In 
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the Lorentz case one can use a similar relation with both a and cr* tending to En i.e. with 
an — 7- En and ct/ — )■ 0. 

Choosing the kernel K as that of the Laplace transform one gets 

<l>(a) = (Q'|e-'^^|g)- J]i?„e-^-. 

n 

This quantity is known to be calculable with the help of the Green function Monte-Carlo 
method at least in the cases when at each total angular momentum and parity values there 
exists not more than one bound state. ^ 

In the framework of no-core shell model calculations other kernels K may probably also 
be used in the present context, see [6l |9]. 



IV. EXAMPLES OF INVERSION 

In this section Eq. ([T]) with Ethr = and with 

^^^^ = ;4[(^lyTiF ^'^^ 

taken as an exact solution is considered. Let us use Eq = 20.7212603615, cf. below. Approx- 
imate solutions denoted fappr will be obtained below. The task here is to find out accuracy 



with which the approximate solutions reproduce the exact one given by Eq. (27). The min- 
imum set of fitting parameters to be used arises when exp['j{E)] is disregarded in Eq. (|6]). 
In the present one-maximum case the corresponding ansatz is 



The parameters C and Ei are expressed in terms of the two other parameters E and /3 with 
the help of the relations (p|) and Then the two parameters, E and (3, are to be fitted. 



Below results pertaining to the approximation (28) are presented 



It is seen that the derivative of the solution (27) lies in the class of the functions (28). 



Therefore if the exact input is used then the two-parameter ansatz (28) should lead to the 



exact solution up to numerical inaccuracies. However, this is not the case when approximate 



At a sufficient accuracy in the calculated {Q'\e '^^\Q) values the quantities i?„ and £"„ can also be 
treated as fitting parameters along with those entering f{E). 
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inputs are employed. Nevertheless, it occurs that the inclusion of additional parameters 
representing the quantity 7(-E) does not lead to noticeable changes of the corresponding 
approximate solutions obtained.^ 

Below some kernels K of interest are considered in Eq. ([T]). The inputs corresponding to 
those kernels are $ = Kf where / is the solution (27). Realistic approximate inputs $appr 
related to those exact $ are generated below and with these $appr approximate solutions 
fappr are subsequently obtained using Eq. Q. 

The average relative deviation of fappr from the exact / is chosen as the criterion of the 
accuracy of a solution. I.e. the quantity 

J_ / f{Ei) — fappr{Ei^ 



X.solution 



i=l 



fiE, 



(IA///I) 



(29) 



is adopted as such a criterion. The sum goes over a large number of E-i values within an 
appropriate range of E. This range is taken to be < < 42 in all the cases, cf. Fig. 1 
below. 



The average relative deviation of ^appr from the exact $ 



X,input 



—y 



<l>(o-i) - ^appr{<yi 



[CTi 



1/2 



(|A<^/<^I) 



(30) 



is taken as the criterion of the accuracy of an approximate input. The sum goes over a large 
number of at values within the ranges of a employed to solve the problem. These ranges are 
specified below. (The quantities $appr(o") have no zeros in the cases considered.) 



At performing the minimization the weight function w{a) entering (10) is chosen to be 



'^/[^apprio-)]'^- I-e. the average square of the relative deviation of K fappr from ^appr is 
minimized. The corresponding fit quality is reported below as well. As the criterion of this 
quality the value 



Xfit 



^ N2 



appr 



{K fappr) icT.i) 



^appr 



1/2 



(l(<^> 



appr 



K fappr) /^apprl) (31) 



j=l \ -appr\ 

is taken. Here fappr is the approximate solution obtained via the minimization procedure of 
Eq. (§. 



^ In some cases changes in the approximate solutions due to inclusion of those parameters prove to be even 
less than numerical uncertainties at finding the minimum. The program "frprmn" from [2, Sec. 10.6, is 
used in the present calculations to search for the minima. 
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A. Lorentz transform 



The normalized Lorentz kernel ( 25 ) 



is used. The parameter aj determines the width of the kernel. At o"/ — the kernel 
turns into 5{(yR — E). The kernel K of the corresponding Eq. ^ to be solved equals 
— TT^^ arctan[(£' — crj?)/cr/]. 

The corresponding input $(o"/j) is given in a range of aR values. If the exact 
is analytic in this range as in the applications of the preceding section then the solution 
corresponding to this $(0"/?) exists, it is unique, and it is independent of the chosen 
range. Indeed, f{E) satisfies the equation that is obtained by the replacement of the lowest 



integration limit Eth^ with —00 if one sets f{E) = at < Eth^. Since the kernel (32) is an- 
alytic f{E) thus defined satisfies the latter equation also in the whole range — 00 < crj? < cxo 
with $(o"r) continued analytically onto this range. And the equation thus obtained is a 
convolution equation having a unique solution. 

Results of the present inversion method are compared below with those emerging from 
the Eq. ^ procedure.^ The basis sets with the E^/"^ threshold behavior entering ^ used 
in almost all practical calculations performed so far are the following^ 

^n{E,a) = E-^-^'^e""^ and a) = E^/V"^/". (33) 

Results that were obtained with these two sets are of similar quality. Below the first of 
the sets is employed. The corresponding calculations will be referred to as the "standard 
inversion" . 

At searching the parameters the Lorentz transforms of the approximate solutions are 
calculated numerically both in the standard inversion and in the present techniques cases. 
The same applies to inversion of the Stieltjes transforms below. 

First let us present the inversion results for the case when the input $(cr) exactly cor- 



responding to the solution (27) is employed. While this input can be calculated directly it 



^ In the Lorentz case the latter procedure was studied in [51 [TSHT7] . In Ref. TB] some other approaches 
were also tested. 



^ More involved basis functions were used in the recent paper 
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is convenient to obtain it in another way aiming to subsequently generate ^appr{o'). Let us 



come to the solution (27) proceeding from the representation of Eq. (17) type. Let the states 



\E'^ in (17) be the eigenstates of the one-particle free-motion Hamiltonian with the orbital 



momentum equal to zero. In this case one may use as the 7 variable in Eq. (17) and the 



summation sign is to be omitted. Let us also set Q = Q' in (17). The quantity f{E) then 



equals |(\E'£;|(5)p. The corresponding wave functions \I'£;(r), normalized to 6{E — E') as in 

1/2 



Eq. IK2\, are 



^i?(r) 



c sin(A;r) 



2M 



(34) 



Here M is the particle mass taken below to be the mass of the nucleon, and E = {likf/{2M). 
Let us also set 



Q(r) 



1 



-•qr 



(35) 



Bound state contributions entering relations of Sec. 3 are to be omitted. This gives 
the expression (27) with Eq = [hr])'^ / {2M) and with the normalization factor replaced by 



4: / (tt EqT]^) . The sum rule value (jSj) and (18) equals 1/(47]^). Below it is set r/ = 1 fm ^. If 
energies are measured in MeV and distances in fm then one gets the Eq value listed above. 



According to the preceding section, see Eq. (24) and below, the transform $ corresponding 



to Eq. (27) can be calculated as 



(r, s)p(ir 



(36) 



where ip is determined from the equation 



2M 



ifj" — sil) = re 



rjr 



(37) 



with 5 = 0"^ + icT/. The solution sought for is localized and satisfies the condition ipi^y s) = 0. 
In all the cases the inputs $ or $appr will be used in the range —2 < ctr < 41.4 with 
aR measured in MeV units. (One may hope to obtain the solution f{E) with a reason- 
able accuracy in the range < < E^ax if the input of a similar accuracy in the range 
—a/ < a < Emax + o'l is employed, cf. also Fig. 1.) 



The exact input deduced from Eqs. (36) and (|37j) is calculated analytically. In Table 1 
the results obtained with this input via the standard inversion are listed at the aj value 
equal to 10 MeV. At this aj value the width of the kernel is not very different from that of 
the solution, cf. Fig. 1 below. Standard inversions were usually done at such a condition. 
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N 


Xfit 


^solution 


5 


8.3-lC^ 


5.2-10-2 


8 


2.3-10-5 


1.4-10-2 


9 


6.2-10-6 


1.3-10-2 


10 


4.1-10-6 


2.65-10-2 



TABLE I: The inversion results. The exact input. The "standard inversion", aj = 10 MeV. 



0-/, MeV 


Xfit 


Xsolution 


2 


1.4-10-"^ 


1.6-10-*^ 


10 


1.45-10-^^ 


2.1-10-*^ 


100 


7.5-10-8 


3.3-10-^^ 



TABLE II: The inversion results. The exact input. The new method, two fitting parameters. 

The results in Table 1 are presented for various choices of the number of basis functions 
retained in the expansion ([2]). 

It is seen that, although the exact input is used, the accuracy of the solution obtained is 
limited and the best accuracy attained is at the level of one per cent. This happens despite 



the fact that the transforms of the basis functions (33) are calculated with high accuracy. 
While the quality of the fit improves monotonically as increases this is not the case as to 
the solution obtained. When N exceeds 9 the quality of the solution deteriorates and the 
instability due to tiny numerical inaccuracies starts to develop. 

In Table 2 the results obtained with the method of the present work in the framework of 



the ansatz (28) are listed. Very different choices of the width parameter ai are considered. 
In all the cases the accuracy of the solution obtained is incomparable with that for the above 
case of the standard inversion. 

Now let us pass to the results at approximate inputs $appr- Let us proceed as follows. 



Eq. (37) is a special case of the dynamic equation {H — a)^! = Q from (23) and (24). Usu- 
ally at solving such equations expansions over basis functions are employed in few-body 
calculations and main inaccuracies in ^appr come from truncations of the expansions. To 



model this situation let us solve Eq. (37) also in the form of a similar expansion. The number 
of terms A^o retained in the expansion will be chosen such that $appr will have the accuracy 
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N 


Xfit 


^solution 


2 


8.0-10-2 


0.22 


3 


2.6-10-2 


0.12 


4 


2.6-10-2 


0.13 


5 


1.35-10-2 


0.71 



TABLE III: The inversion results. An approximate input. The "standard inversion", aj = 10 MeV. 



at a per cent level. Specifically, let us use the approximation 



No 



^{r, ~ 5^ Cn{s)Mr), Mr) = [n{n + l)]-'/%-'/W^_,{r/b) exp[-r/(26)]. (38) 



n=l 



Here b = 0.3 fm-^ is the length parameter, and Lf^{x) are the Laguerre polynomials. The 
basis functions are orthonormalized. The expansion coefficients c„ are determined from the 
projection linear equations. The corresponding matrix elements entering the equations and 



the resulting approximate inputs (36) are calculated analytically. The same is done below 



in the case of the Stieltjes transform. 



First let us consider again the aj = 10 MeV case. Let us set A^o = 10 in Eq. (38). This 



leads to the corresponding ^appr that has a three per cent accuracy. The accuracy is defined 



as in Eq. (30 ) 



Xinput = 3.0 ■ 10-2. (39) 

The results of the standard inversion are presented in Table 3. One sees that at a three 
per cent accuracy of the input the best accuracy of the solution is twelve per cent and thus 
is unsatisfactory. (One should mention at the same time that many good accuracy results 
were previously obtained in such type problems applying the standard inversion. Probably 
the accuracy of the input was considerably higher than three per cent in those calculations.) 
Furthermore, it occurs that the approximate solutions at Nmax = 3 and Nmax = 4 are very 
close to each other. However, they are still not close to the true solution. Thus, stability 
with respect to the number of basis functions retained at performing the standard inversion 
is not sufficient for an approximate solution to be close to the true solution. Stability with 
respect to variations of ^appr is also required. 

Now let us pass to the result obtained with the method of the present paper in the same 
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case (see Eq. (39) and above it, cx/ = 10 MeV). One has 

Xfit = 2.9 ■ 10-^ Xsoiuuon = 1.65 ■ 10-^. (40) 



The accuracy of the solution obtained is even better than the accuracy (39) of the input in 
contrast to the standard inversion resuhs. 

Perhaps the latter results might be somewhat improved due to use of basis functions 
decreasing as powers of E at large E values. This would require inclusion more non-linear 
parameters in each of basis functions. Anyway, the advantage of the present method is that 
in its frames there is no need to worry about the choice of the regularization parameter and 
instability does not arise at all. 

Let us also list the results obtained with the present method at other o"/ values. Let 
us first consider the case of a very narrow kernel, o"/ = 2 MeV. In this case the decrease 
of iplr) as r increases is slow with many oscillations present in its range of localization. 



Therefore many basis functions in (38) are required to reproduce ijj{r) with a reasonable 
accuracy. With A^o = 60 in (38) one gets Xmput = 2.6 • 10"^. Using the corresponding ^appr 



one obtains with the present method 

Xfu = 2.6 ■ 10-2, Xsoiution = 8.6 ■ 10-^ (41) 

The accuracy of the solution obtained is even considerably higher than that of the input in 
this case. 

Let us now consider the opposite case of a very broad kernel, aj = 100 MeV. Very few 



basis functions in (38) are to be retained in this case to provide a reasonable accuracy of 



the input. At the number A'o = 3 of these basis functions one gets 

Xinput = 1.55 • lO-\ (42) 

The standard inversion is absolutely inapplicable in this case since the instability develops 
even for the lowest A^ values in (|2]). However, the method of the present work allows to cope 
with this case as well. Using the corresponding ^appr one obtains with this method 

Xfit = 4.8 ■ 10-^ Xsoiution = 7.7 ■ lO-l (43) 
The attained accuracy of the solution is still very reasonable. 
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The results for this case are displayed in Fig. 1. The dashed curve represents the exact 
solution and the dotted curve represents the approximate one. They are practically indis- 
tinguishable. The solid curve represents the approximate input employed at performing the 
inversion. It is practically a constant. However, its tiny deviations from a constant prove to 
allow reconstructing the solution with a very good accuracy. 



a, = 100 MeV 



E or a (MeV) 



FIG. 1: The inversion results. The new method, two fitting parameters. Accuracy of the approx- 



imate input is given by Eq. (42). This input is represented by the solid curve. The dashed curve 



and the dotted curve represent, respectively, the exact solution and the approximate one. 



If at aj = 100 MeV the number A^o of the basis functions in (38) is increased up to e.g. 
A^o = 10 one gets 

Xinput = 9.3 ■ 10-^ (44) 



Using the corresponding ^appr one obtains in this last case with the present method 



Xfu = 2.1 -10- 



Xsolution 2.1 ■ 10 



(45) 



This is to be compared with the values (40) obtained at the same Nq. One sees that, at the 



same number of basis functions employed to solve the dynamics equation, use of a very broad 
kernel is more preferable in the problem than use of a kernel with the width comparable 
with that of the solution. This is due to faster convergence with Nq at solving the equation 



(37) of the type of {H — a)"^ = Q, see the values in Eqs. (44) and (39) 



Comparison between results obtained from calculations at different aj may serve for 
checking purposes. 
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B. Stieltjes transform 



The Stieltjes kernel is 

K{s,E) = {E - s)-\ (46) 

Eq. ([T| is solved with given in a range of real negative s values. If the exact is 
analytic in this range as in the applications of Sec. 3 then there exists a unique solution 
corresponding to this that is independent of this s range. This follows from the fact that 
the Stieltjes transform is the iterated Laplace transform. Eq. ^ with the kernel — ln{E — s) 
is solved. 



The same procedure as above that involves Eqs. (34) and (35) is used to calculate the 



exact input and approximate ones. As it is seen from the preceding section, Eq. (23), the 



exact input $ that corresponds to the solution (27) can be calculated in the Stieltjes case 
as 



= / ip{r,s)re~'^''dr (47) 
Jo 

where '4'{r,s) is the localized solution to Eq. ([37)) at the condition ip{0,s) = 0. The same 



units and the same rj value as above are used. The same procedure of expanding over the 



basis set (38) is used to obtain approximate inputs. 

At performing the inversion the inputs in the ranges Smax > s > Smin are employed with 
various Smax and with Smin = Smax — 41.4 MeV. To this aim, it is expedient that |s| is 
confined to values not too much exceeding those E values at which f{E) is substantially 
different from zero, cf. Fig. 1. At larger |s| values the transform behaves as const ■ and 
it contains little additional information on f{E). As to the choice of the Smax value, the 
transform exists for s < only and convergence of the approximate transform deteriorates 
when s approaches zero. 

At Smax = —2 MeV in the case of the exact input one obtains 

X/it = 3.45 ■ 10-^2^ Xsoiution = ^-9 ■ lO-'\ (48) 



In the case of an approximate input obtained with Nq = 5 in (38) one has Xmput = 1-3 ■ 10 ^. 
Performing inversion with the corresponding ^appr one gets 

X/it = 5.5- 10-^ Xsoiution = 0-10- (49) 
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While this result is unsatisfactory, the increase of the number A^o of the basis functions up 
to Nq = 7 leads to 

Xinput = 2.2 ■ 10~\ Xfit = lA-10-\ Xsolution = 1-5 -10-^. (50) 

Let us mention that for obtaining the solution with a similar accuracy when the standard 
inversion is used the accuracy of the input should be at the level of 10^^-10"^. 

If at the same s^ax = ~2 MeV value one increases further the number of retained basis 
functions up to A^o = 10 one obtains 

Xinpni = 1.2 ■ 10-^ X/ii = 1-4 ■ 10-^ Xsolution = l.l ■ l^-^ . (51) 

If at the same A^o = 10 value one uses Smax = — 10 MeV and Smax = ~20 MeV one obtains, 
respectively, 

Smax = -10 MeV, Xinput = 1-4 " 10"^ Xfit = 9.9 ■ 10-^ Xsolution = 2.2 ■ 10-^ (52) 

Smax = -20 MeV, Xinput = 2.3 • 10-^ Xfit = 1.5 ■ 10-^ Xsolution = LO • 10"'. (53) 

Thus, despite the fact that the ratio Xsolution/ Xinput becomes less favorable as Smax decreases, 
the accuracy of the solution itself increases. The reason for this is similar to that in the 
Lorentz transform case. 

Comparison of the results obtained from calculations at different Smax values may serve 
for checking purposes. Besides, use of the Lorentz and Stieltjes transforms simultaneously 



for the same purposes is convenient. Indeed, the dynamic equations of Eq. (37) type are the 
same in both cases. In addition, in the Lorentz case it is in general convenient to calculate 



$ from the Stieltjes type relations, see Eq. (26). 



C. Laplace transform 

When conventional regularization methods are used to invert the Laplace transform the 
instability problem is known in general to be severe. The Laplace transform considered here 
is 

/■oo 

/ e-^^/(E)rfE = $(^), /(O) = /(oo) = 0. 
Jo 

Inversion is performed in the form 

POO 

/ e-'^f'{E)dE = z<^{z). 
Jo 
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The right-hand side values in the range < z < Zmax with Zmax — 1-9304 MeV~^ are used 
when performing the inversion. If is a typical scale pertaining to the solution f{E), cf. 
Fig. 1, then at z ^ 1/E the transform takes its asymptotic form and contains little additional 
information on f{E). In the present case only minimization on grids is performed without 
a subsequent conventional minimization. 



The transform ^{z) of the exact f{E) of Eq. (27) is calculated numerically. The inversion 



with this accurate ^{z) gives (using the two-parameter ansatz (28)) 

Xf^t = 2.6 ■ 10-^ Xsolution = 8.0 ■ 10^^ 

In the literature, in many-body calculations inputs to Laplace transforms are obtained 
with the help of Monte-Carlo integrations. Therefore it will be assumed here that inaccu- 
racies of an input are of a random nature. Random inaccuracies may be modeled in various 
ways. Below it is done as follows. At each z value the quantity v^{z) is added to the 
accurate input ^{z) where t; is a random variable such that its mean value is zero. Let us 
denote its dispersion, i.e. {v"^)- One may set v = tq where q is a, random variable with 
the unit dispersion. It is assumed that the distribution of q is the normal distribution and 
the values of q are taken randomly in accordance with this distribution [7J at each z = Zi 
value. 

The average relative error r of the input is taken to be 5.0 ■ 10~^. Inversion with the 
approximate input ^appr{,z) thus generated gives 

Xf^t = 4.9 ■ 10-2, Xsolution = 8.9 ■ 10"^ 

The value of Xfu is such as one could expect. The result for XsoiuUon obtained may be 
considered definitely good taking into account rather large relative error in the input. 



V. CONVERGENCE OF THE METHOD 



In this section it is shown that if both the quantity ||$ — $appr || is small and saturation 
is achieved with respect to the number M of parameters retained e.g. in Eq. ([T]) then the 
approximate solution /m from Eq. ^ is necessarily close to the true / everywhere except 
perhaps for the points of maxima and minima of /. The case when the uncertainty $ — ^appr 
is not random is considered here so that one can speak of the norm ||$ — $appr ||- (In the 
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above example the method works well also in the case of random uncertainty.) Besides, the 
influence of round-off errors on the minimization procedure ^ is disregarded below. 

1. First, let us show the following. Let /m be the solution to Eq. If ||$ — $appr. || 
is small and saturation with respect to the number M of parameters determining Jm is 
achieved then the quantity ||$ — K/mW = \\K{f — /m)|| is small as well. Namely, 

<2||<f -$„pp,||+ei (54) 

where ei may be done arbitrarily small due to increase in the number M of parameters 
retained. 

Indeed, one has 

||^(/-/m)|| < - ^apprW + W^appr - K/mW- (55) 

Let f^j be an arbitrary function of the structure of Eqs. ^ and ([t]) with the same M 
parameters retained as in the case of /m- Due to Eq. Q one can write 

W^appr - KJmW < ||$appr - ^/a/II < 1 1 "^appr " $ 1 1 + 1 1 $ " ^/m 1 1 • (56) 

Let us choose f^j to be arbitrarily close to exact / taking M value to be sufficiently large. 
Since the K operator is continuous, at any ci > one can find such 6 in the corresponding 



relation ||/ — < S that one gets ||$ — Kf^j\\ < ei. Then combining Eqs. (55) and (56) 
one obtains Eq. (54). (The above norm in the / space may be deflned differently from that 
in the $ space.) 

2. Now the following question is to be addressed. Suppose that \\K{f — /a/)]] is small. 
What are implications of this fact for the behavior of / — /m itself? 
Below let us deflne the K operator as follows, 

{K^)ix)= [ K{x,y)^Piy)dy, c<x<d (57) 

J a 

where Ktp G i^2(c, d), and the functions ip belong to the class of piecewise continuous func- 
tions. The notation ||F|| will refer to the L2 norms below. Performing a change of variables 
the norm of Eq. ( 10 ) may be rewritten as the L2 norm. Let us also assume that K is such 
that if ll-R'^/'ll = then ||-?/^|| = in the class of piecewise continuous ip. (This condition is 
equivalent to the condition that the solution to the equation Kf = $ is unique also in the 
class of piecewise continuous /.) 
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Consider the mapping Kz = u, where z{y) are smooth uniformly bounded functions, 
||z|| < Z, where ||z|| is the L2{a,b) norms of z{y). Let be the L2{c,d) norm of u{x). 

The following will be proved. At any yi and y2 belonging to [a, b] and at any e > one 
can find such S that the relation 



1/2 



z{y)dy 



yi 



< e (58) 



holds true for all z satisfying the condition < 6. (Here u = u[z].) One can split the 



a, 6] integration domain from Eq. (57) into segments of a given length A = y2 — yi- The 



above mentioned 6 value may be found such that Eq. (58) fulfills simultaneously for all these 
segments. 



Eq. (58) provides a tool for establishing the local properties of the approximate solution. 
Consequences for our purposes of this relation will be discussed below after its derivation. 
The above condition that z = f — /m are uniformly bounded is equivalent to the condi- 
tion that /m are uniformly bounded. If required, the latter condition may be imposed at 
performing the minimization procedure. 



To obtain the relation (58) let us write 



"b rb rd 

|2 



u\\'= / z{y)}C{y,y')z{y')dydy' with IC{y,y') = K{x,y)K{x,y')dx, (59) 



O' < y,y' < b. In addition it will be assumed here that | |/C| p = J^^[/C(y, y')Ydydy' is finite. 
When a or 6 is infinite this may be not the case. Then for the purpose of the present 
reasoning the corresponding integration limit may be replaced with a large finite number. 
Indeed, such a replacement does not influence physics, and the magnitude of the norm ||/C|| 
will be of no relevance below. 

For this purpose in the above considered cases of the Lorentz, Stieltjes and Laplace 
transforms one would need to replace the infinite upper integration limit with a finite R 
value. It can be seen that in those cases such a replacement leads merely to the replacement 
of the solution fiy) with f{y)9{R — y). It can be seen that the above mentioned property 
I I = ^ \ \ = 0, being piecewise continuous, is valid in the case of those transforms 
both in their original form and after the replacement of the upper integration limit with R. 

® Below to employ expansions over continuum spectrum eigenfunctions pertaining to these transforms could 
be an alternative. 
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The kernel /C is the Fredholm one and is symmetric. Therefore it possesses eigenfunctions 
and it is possible to number them [12]. Let us denote them (pn, n = 1, . . ., and let us denote 
/x„ the eigenvalues, /C0„ = /in0n- Below let us choose 0„ to be orthonormalized. 

Note that zero cannot be an eigenvalue of /C, and /i„ > at any n. Indeed, according to 

■ \\K(j)n\\'^, and the kernel K is such that ||0„|| 7^ ^ ll-^0n|| 7^ 0. 



Eq. (59) /i„ 



The set of the eigenfunctions of the kernel K, is complete in the class of piecewise contin- 
uous functions in the sense of the approximation in the L2 norm. 

Indeed, for any function ip orthogonal to all 0^ the equality Kip = is valid almost every- 
where in [a,b] [IH]. Hence {ip,]Cip) = for such ip where (x, v^) denotes the standard scalar 



product. But according to the definition (59) the relation {ip,}Cip) = yields = 

which was adopted above to be not possible at \\ip\\ 7^ in case of the kernels K under 
consideration. Thus if ip is orthogonal to all 0„ then \\ip\\ = 0. I.e. the set {(pn} is closed 
and therefore it is complete. The set is thus infinite. 



Let us pass now now to Eq. (58). Let c„ be the Fourier coefficients in the expansion of z 
over (pn, 

00 

z{y) --^Cn(pn{y)- (60) 



n=l 



Let us represent the integral in Eq. ( 58 ) as 



/ z{y)dy = ^Cndn (61) 

•^2/1 n=l 

where dn are the Fourier coefficients in the expansion of 9 {1/2 ^ y) — Gijji — y) over (pn{y)- 
Eq. (61) follows from the completeness property of the set {(pn}- Since all ||2;|| do not exceed 



some Z value one has at any no 



V c 



< z 



_ra=?iO+l 



1/2 



(62) 



n=no+l 

Taking this into account along with the fact that the function 6 {1/2 — y) — 6{yi — y) has a 
finite norm one sees that at any e > one can find such uq that the relation 



V c 



< 



(63) 



n=no+l 

is valid. The estimate is uniform with respect to z{y) out of the class considered. 
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Furthermore, one has 



n=l 



<Si<S2<S3 



(64) 



where 



^1 



no 
n=l 



1 1/2 



5. 



^3 = 



no 



(Z/2 - 1/1 
1/2 



,1/2 



E2 2 



n=l 



(1/2 -Z/l)^^ 



|m||(?/2 - 2/1 



-7„^l/2 



and fir. 



[no) 



min(/ii,/X2, 



f^min (^o) 

/x„q). Let US set 

(^o) 



u < 



2(y2-yi)V2- 



(65) 



Combining then Eqs. (63) and (64) one obtains that the absolute value of the quantity (61) 



does not exceed e, i.e. at arbitrary e one comes to Eq. (58). 



If the [a,b] integration domain in Eq. (57) is split into the segments of the length of 
A = y2 — yi then at a given e one can adopt the highest of the riQ values pertaining to these 



segments as the common uq. The value in the estimate (65) can be chosen the same 



for all the segments. Then Eq. (58) will be valid for all the segments simultaneously. 



3. The property (58) shows that at small the quantity z may be not small only 



because of narrow peaks of high amplitude. And if the requirement is imposed that the 
approximate solution Jm has the same number of maxima and minima as the true / has 
such narrow peaks of high amplitude are forbidden everywhere except for the points of 
maxima and minima of /. Therefore z is necessarily small everywhere except possibly for 
these points. 

Let us show this in more detail. Suppose that the [a, b] integration domain is split into 
segments of the length of A. Below A is considered to be sufficiently small. Suppose that 



the property (58) is fulfilled for all the segments simultaneously. This imposes the following 



limitation on the behavior of z{y). It varies within the band \z\ < e/A and out of the 
band peaks with widths not exceeding 2A are only permissible. (I.e. if a segment [yi,y2] 
exists at which z{y) > e/A or z{y) < — e/A for all the y values belonging the segment then 



necessarily y2 — yi < 2 A.) Indeed, otherwise the condition (58) would be violated at some 
of the segments. 
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Let us exclude ±A vicinities of maxima and minima of / from the consideration (cf. 
the beginning of Sec. 2 in this connection). Let us impose now the requirement that the 
approximate solution Jm has the same number of maxima and minima as the true / has. 
This requirement puts limits on the amplitudes of the above mentioned peaks possibly 
present beyond the band. These amplitudes cannot exceed the |/'| ■ 2 A + e/A values where 
/' are average derivatives of the true solution in those segments (of the length not exceeding 
2A) where the peaks are located. Indeed, otherwise, due to these peaks, the quantity 
I'm = f — z would include more maxima and minima than / does. But in this consideration 
A may be chosen arbitrarily small from the beginning. And when A is given e(A) may be 
taken arbitrarily small provided that \\K{f — /m)|| is sufficiently small. At this condition 
\z\ becomes smaller than any prescribed value. 

Thus in the whole [a, b] integration domain with the points of maxima and minima of 
/ excluded the following is valid. If the mentioned requirement on the number of maxima 
and minima is imposed then there exist 6i values such that when — /m)|| < ^i, cf. 



Eq. (54), the max|/ — /m] quantity is necessarily smaller than any prescribed value. 

In conclusion, techniques to solve integral equations of the first kind with an approximate 
input are proposed. It is proved that, at the conditions and reservations listed above, the 
present method provides the solution stable with respect to perturbations of an input. No 
regularization is required at solving a problem in this way. The fact that one need not deal 
with a regularization parameter may in practice lead to a higher accuracy. Some inversion 
problems allowing comparison with an exact solution of a rather simple structure are studied 
numerically. They include inversions of the Lorentz, Stieltjes and Laplace transforms and 
involve systematic and random errors of an input. The results prove to be very satisfactory. 
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